</div>

Tarea 3

Física Computacional

</div>

Autor: Celeste Castro Granados $\mathbb C \hat{e}l \mathbb s$

Primero llamaremos a las bibliotecas que vamos a necesitar:

In [1]:
using Plots
using LaTeXStrings

1. En clase vimos el método de Euler mejorado en el que nosotros introducimos el Jacobiano de forma analítica para ser evaluado por el método de integración. Realice una función que implemente el mismo método de integración, pero en lugar de calcular el Jacobiano con una función que lo evalúa analíticamente, ahora sea calculado de manera numérica utilizando la derivación numérica.

En general, una ecuación diferencial se representa por:

$$ \frac{d\vec x}{dt} = \vec g (\vec x, t) $$

Y sabemos por lo visto en clase que el método de Euler mejorado realiza una aproximación mediante la derivada del sistema de ecuaciones utilizando la matriz jacobiana y el desarrollo de Taylor de la $\vec x (t)$ por medio de la implementación del siguiente algoritmo:

$$ \vec x_{n+1} = \vec x_n + \delta \vec g ( \vec x_n , t_n ) + \frac{1}{2} \delta^2 \mathbb{J}^* ( \vec x_n, t_n) \cdot \vec g^* ( \vec x_n, t_n ) $$

En donde: $\mathbb{J}^* = \big ( \mathbb{J}, \frac{\partial \vec g}{\partial t} \big )$ con $\mathbb{J}$ el Jacobiano, y $\vec g^* = (\vec g, 1)$. (se agrega un 1 como entrada extra en $g^*$ para que no haya problemas al momento de multiplicar por $\mathbb{J}^*$, i.e. que las dimensiones concuerden)

Ahora, queremos una función que calcule de manera numérica la matriz $\mathbb{J}^*$. Para esto, recordamos que:

$$\mathbb{J}^* = \begin{bmatrix} \frac{\partial g_1}{\partial x_1} & \frac{\partial g_1}{\partial x_2} & ... & \frac{\partial g_1}{\partial x_n} & \frac{\partial g_1}{\partial t} \\ \frac{\partial g_2}{\partial x_1} & \frac{\partial g_2}{\partial x_2} & ... & \frac{\partial g_2}{\partial x_n} & \frac{\partial g_2}{\partial t} \\ . & & & \\ . & & & \\ . & & & \\ \frac{\partial g_n}{\partial x_1} & \frac{\partial g_n}{\partial x_2} & ... & \frac{\partial g_n}{\partial x_n} & \frac{\partial g_n}{\partial t} \end{bmatrix}$$

Por tanto, necesitamos calcular las derivadas parciales de la función $\vec g$. Dichas parciales las calcularemos mediante la siguiente expresión:

$$ \partial_{x_j} g_{i,n} = \bigg (\frac{1}{2h} \bigg ) \bigg [ g_i \big (( x_{1,n}, ..., x_{n,n}, t_n ) + (0,...,h,...,0) \big ) - g_i \big (( x_{1,n}, ..., x_{n,n}, t_n ) - (0,...,h,...,0) \big ) \bigg] $$

En donde, $h$ se encuentra en la j-ésima posición ya que estamos calculando la parcial respecto a $x_j$.

A continuación, procedemos a programar la función que nos de esta matriz:

In [2]:
function jacobiano_t(xi,h,f,t)
    #xi-consiste en las variables del problema, donde calcularemos las parciales
    #h-paso utilizado para obtener las derivadas parciales
    #f-función a la cual sacaremos el jacobiano
    #t-tiempo
    
    #Definimos la matriz que vamos a ir llenando:
    jacob=zeros(length(f(xi,t)),length(xi)+1)
    
    for i in 1:length(f(xi,t)), j in 1:length(xi) #ciclo para ir llenando las entradas de la matriz
        #xp1,xp2-vectores que tienen un aumento o disminución en la i-ésima entrda.
        xp1=zeros(length(xi))
        xp2=zeros(length(xi))
        
        for k in 1:length(xi) #ciclo con el cual definiremos las entradas de los vectores donde evaluaremos la función para obtener la derivada parcial
            if  k==j #incremento a la j-ésima entrada del vector
                xp1[k]=xi[j]+h
                xp2[k]=xi[j]-h
            else
                xp1[k]=xi[j]
                xp2[k]=xi[j]
            end
        end
        
        jacob[i,j]=(f(xp1,t)[i]-f(xp2,t)[i])/(xp1[j]-xp2[j]) #entrada ij de la matriz, sin tomar en cuenta la entrada de la parcial temporal
        jacob[i,length(xi)+1]=(f(xi,t+h)[i]-f(xi,t-h)[i])/(2h) #componentes de la columna de la parcial temporal
    end
    return jacob
end
Out[2]:
jacobiano_t (generic function with 1 method)

Luego, como ya definimos la función que calcula de manera numérica la matriz $\mathbb{J}^*$, modificamos el algoritmo de Euler mejorado que hace uso del jacobiano con dependencia temporal visto en clase para integrarle nuestra función:

In [3]:
function Euler_mejorado_Jacobnum(edo,edo1,p_ini,t)
    #edo-función de la ecuacion diferencial
    #edo1-función de la ecuacion diferencial con un elemento 1 en la última entrada
    #p_ini-condiciones iniciales de la posición
    #t-tiempo
    # J y g con dependencia temporal.
    sol = zeros(length(t),length(p_ini))
    sol[1,:] = p_ini
    δ = t[2]-t[1]
    for i in 1:(length(t)-1)
        eval_edo = edo(sol[i,:],t[i])
        g=edo1(sol[i,:],t[i]) #función por la que multiplicaremos el jacobiano
        if length(eval_edo) == length(p_ini)
            sol[i+1,:] .= sol[i,:] .+ δ .*eval_edo
        else
            sol[i+1,:] .= sol[i,:] .+ δ .*eval_edo[1:(end-1)]
        end
        sol[i+1,:] .+= 0.5*(δ^2) .*(jacobiano_t(sol[i,:],δ,edo,t[i]) *g )
    end
    return sol
end
Out[3]:
Euler_mejorado_Jacobnum (generic function with 1 method)

2. Utilizando los integradores de Euler mejorado con Jacobiano analítico, Euler mejorado con Jacobiano numérico, Runge-Kutta de $2^°$ orden y Runge-Kutta de $4^°$ orden, encuentre un paso de integración que proporcione una precisión de $10^{-4}$ en la solución para la siguiente ecuación:

$$\dddot{y} + \ddot{y}^2 - 3\dot{y}^3 + cos^2 y = e^{-t} sin (3t)$$

Con condiciones iniciales $\ddot{y}(1)=1$, $\dot{y}(1)=2$, $y(1)=1$. En particular que sucede en el intervalo $t \in [1,2]$ y $t \in [1,2.1]$.

Recordamos primero los algoritmos vistos en clase de los métodos que utilizaremos:

Euler mejorado con Jacobiano analítico

In [4]:
function Euler_mejorado_J_t_p(edo,p_ini,t,jacob)
    # J y g con dependencia temporal
    sol = zeros(length(t),length(p_ini))
    sol[1,:] = p_ini
    δ = t[2]-t[1]
    for i in 1:(length(t)-1)
        eval_edo = edo(sol[i,:],t[i])
        if length(eval_edo) == length(p_ini)
            sol[i+1,:] .= sol[i,:] .+ δ .*eval_edo
        else
            sol[i+1,:] .= sol[i,:] .+ δ .*eval_edo[1:(end-1)]
        end
        sol[i+1,:] .+= 0.5*(δ^2) .*(jacob(sol[i,:],t[i]) *eval_edo )
    end
    return sol
end
Out[4]:
Euler_mejorado_J_t_p (generic function with 1 method)

Runge-Kutta de $2^°$ orden

In [5]:
function RK_2(edo,x_ini,t)
    sol = zeros( length(t) , length(x_ini) )
    sol[1,:] .= x_ini
    δ = t[2]-t[1]
    for i in 1:length(t)-1
        k1 = sol[i,:] .+ 0.5*δ .*edo(sol[i,:],t[i])
        sol[i+1,:] .= sol[i,:] .+ δ*edo(k1,t[i]+0.5*δ)
    end
    return sol
end
Out[5]:
RK_2 (generic function with 1 method)

Runge-Kutta de $4^°$ orden

In [6]:
function RK_4(edo,x_ini,t)
    sol = zeros( length(t) , length(x_ini) )
    sol[1,:] .= x_ini
    δ = t[2]-t[1]
    for i in 1:length(t)-1
        k1 = edo(sol[i,:],t[i])
        k2 = edo(sol[i,:] .+ 0.5*δ.*k1 , t[i] + 0.5*δ)
        k3 = edo(sol[i,:] .+ 0.5*δ.*k2 , t[i] + 0.5*δ)
        k4 = edo(sol[i,:] .+ δ.*k2 , t[i] + δ)
        sol[i+1,:] .= sol[i,:] .+ (δ/6.0).*(k1 .+ 2.0.*k2 .+ 2.0.*k3 .+ k4)
    end
    return sol
end
Out[6]:
RK_4 (generic function with 1 method)

Euler mejorado con Jacobiano numérico

En este caso, utilizaremos el algoritmo desarrollado en el problema 1.

Definición de EDO y Jacobiano analítico

Vamos a analizar ahora la ecuación que tenemos:

Para facilitar las cosas, podemos ver a $y$ como la posición:

$$ \Rightarrow \dot{y}= v \quad (\text{velocidad}) $$$$ \Rightarrow \ddot{y}=\dot{v}= a \quad (\text{aceleración}) $$

Por lo que, la ecuación diferencial se transforma en:

$$ \dot{a} + a^2 -3v^3 + cos^2 y = e^{-t} sin(3t) $$

Y tenemos las siguientes 3 EDO's de 1er orden:

$$ \dot{a} = e^{-t} sin(3t) - a^2 +3v^3 - cos^2 y $$$$ \dot{v} = a $$$$ \dot{y} = v $$

Y utilizando las condiciones iniciales dadas sabemos que: $a(1) =1$ , $v(1)=2$, $y(1)=1$.

Por lo tanto, podemos definir las funciones g's como:

$$ g_1 = e^{-t} sin(3t) - a^2 +3v^3 - cos^2 y $$$$ g_2 = a $$$$ g_3 = v $$

Luego, calculamos de manera analítica el Jacobiano:

$$\mathbb{J}^* = \begin{bmatrix} \frac{\partial g_1}{\partial a} & \frac{\partial g_1}{\partial v} & \frac{\partial g_1}{\partial y} & \frac{\partial g_1}{\partial t} \\ \frac{\partial g_2}{\partial a} & \frac{\partial g_2}{\partial v} & \frac{\partial g_2}{\partial y} & \frac{\partial g_2}{\partial t} \\ \frac{\partial g_3}{\partial a} & \frac{\partial g_3}{\partial v} & \frac{\partial g_3}{\partial y} & \frac{\partial g_3}{\partial t} \end{bmatrix}$$$$ \Rightarrow \mathbb{J}^* = \begin{bmatrix} -2a & 9v^2 & -2 sin(y) cos(y) & e^{-t} (2 cos (3t) -sin (3t)) \\ 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \end{bmatrix}$$

Por lo tanto, ya tenemos todo listo. A continuación, procedemos a programar la ecuación diferencial y el Jacobiano:

In [7]:
function edo_y_1(variables,t) 
    #variables=(a,v,y)
    g1=(exp(-t)*sin(3*t))-(variables[1]^2)+3*(variables[2]^3)-((cos(variables[3]))^2)
    g2=variables[1]
    g3=variables[2]
    return [g1,g2,g3,1] #Se agrega un 1 para que las dimensiones coincidan y se pueda multiplicar por el jacobiano
end
Out[7]:
edo_y_1 (generic function with 1 method)
In [8]:
function jacob_edo_y(variables,t)
    Jacob = zeros(3,4) #definimos una matriz de ceros de 3x4 para poder llenarla
    Jacob[1,1] = -2*variables[1]
    Jacob[1,2] = 9*variables[2]^2
    Jacob[1,3] = -2*sin(variables[3])*cos(variables[3])
    Jacob[1,4] = exp(-t)*(2*cos(3*t)-sin(3*t))
    Jacob[2,1] = 1
    Jacob[3,2] = 1
    return Jacob
end
Out[8]:
jacob_edo_y (generic function with 1 method)

Y recordamos que también debemos definir la función de la ecuación diferencial sin el 1 extra en el arreglo que entrega como resultado:

In [9]:
function edo_y(variables,t) 
    #variables=(a,v,y)
    g1=(exp(-t)*sin(3*t))-(variables[1]^2)+3*(variables[2]^3)-((cos(variables[3]))^2)
    g2=variables[1]
    g3=variables[2]
    return [g1,g2,g3]
end
Out[9]:
edo_y (generic function with 1 method)

Finalmente, definimos las condiciones iniciales para poder utilizarlas al momento de calcular las soluciones con cada uno de los métodos:

In [10]:
cond_inic_y=[1,2,1]; #condiciones iniciales

A continuación, vamos a obtener las soluciones con cada uno de los métodos mencionados arriba, encontrando en cada caso el paso adecuado para el cual obtenemos una precisión de $10^{-4}$, i.e. que obtenemos al menos 4 cifras significativas:

Euler con Jacobiano analítico

Vamos a ver que sucede primero en el intervalo [1,1.8]:

In [11]:
t_euler_a1= collect(1:0.001:1.8)
solucion_euler_a1 = Euler_mejorado_J_t_p(edo_y_1,cond_inic_y,t_euler_a1,jacob_edo_y)
Out[11]:
801×3 Array{Float64,2}:
   1.0       2.0      1.0    
   1.02275   2.00101  1.002  
   1.0455    2.00205  1.004  
   1.06823   2.0031   1.006  
   1.09096   2.00418  1.00801
   1.11368   2.00528  1.01001
   1.13639   2.00641  1.01202
   1.15908   2.00756  1.01403
   1.18177   2.00873  1.01603
   1.20445   2.00992  1.01804
   1.22712   2.01114  1.02005
   1.24977   2.01237  1.02207
   1.27242   2.01364  1.02408
   ⋮                         
 382.429    37.0242   6.99894
 388.496    37.4096   7.03616
 394.693    37.8012   7.07376
 401.022    38.199    7.11176
 407.487    38.6032   7.15016
 414.091    39.014    7.18897
 420.839    39.4314   7.22819
 427.735    39.8557   7.26783
 434.782    40.2869   7.3079 
 441.985    40.7252   7.3484 
 449.348    41.1709   7.38935
 456.875    41.6239   7.43075

Notamos que en el intervalo [1,1.8], las soluciones se encuentran bien definidas, sin embargo, veremos que si nos movemos a los intervalos [1,1.9], [1,2] ó [1,2.1] las soluciones ya no se encontrarán bien definidas:

In [12]:
t_euler_a2= collect(1:0.001:2.1)
solucion_euler_a2= Euler_mejorado_J_t_p(edo_y_1,cond_inic_y,t_euler_a2,jacob_edo_y)
Out[12]:
1101×3 Array{Float64,2}:
   1.0        2.0        1.0    
   1.02275    2.00101    1.002  
   1.0455     2.00205    1.004  
   1.06823    2.0031     1.006  
   1.09096    2.00418    1.00801
   1.11368    2.00528    1.01001
   1.13639    2.00641    1.01202
   1.15908    2.00756    1.01403
   1.18177    2.00873    1.01603
   1.20445    2.00992    1.01804
   1.22712    2.01114    1.02005
   1.24977    2.01237    1.02207
   1.27242    2.01364    1.02408
   ⋮                            
 NaN        NaN        NaN      
 NaN        NaN        NaN      
 NaN        NaN        NaN      
 NaN        NaN        NaN      
 NaN        NaN        NaN      
 NaN        NaN        NaN      
 NaN        NaN        NaN      
 NaN        NaN        NaN      
 NaN        NaN        NaN      
 NaN        NaN        NaN      
 NaN        NaN        NaN      
 NaN        NaN        NaN      

Euler con Jacobiano numérico

Nuevamente veamos primero que sucede en el intervalo [1,1.8]:

In [13]:
t_euler_n1 = collect(1:0.00001:1.8)
solucion_euler_n1= Euler_mejorado_Jacobnum(edo_y,edo_y_1,cond_inic_y,t_euler_n1)
Out[13]:
80001×3 Array{Float64,2}:
   1.0       2.0      1.0    
   1.00023   2.00001  1.00002
   1.00046   2.00002  1.00004
   1.00068   2.00003  1.00006
   1.00091   2.00004  1.00008
   1.00114   2.00005  1.0001 
   1.00137   2.00006  1.00012
   1.00159   2.00007  1.00014
   1.00182   2.00008  1.00016
   1.00205   2.00009  1.00018
   1.00228   2.0001   1.0002 
   1.0025    2.00011  1.00022
   1.00273   2.00012  1.00024
   ⋮                         
 456.131    41.5793   7.42659
 456.207    41.5838   7.42701
 456.283    41.5884   7.42742
 456.359    41.593    7.42784
 456.435    41.5975   7.42825
 456.511    41.6021   7.42867
 456.587    41.6066   7.42908
 456.663    41.6112   7.4295 
 456.74     41.6158   7.42992
 456.816    41.6203   7.43033
 456.892    41.6249   7.43075
 456.968    41.6295   7.43117

Y ahora para el intervalo [1,2.1]:

In [14]:
t_euler_n2= collect(1:0.00001:2.1)
solucion_euler_n2= Euler_mejorado_Jacobnum(edo_y,edo_y_1,cond_inic_y,t_euler_n2)
Out[14]:
110001×3 Array{Float64,2}:
   1.0        2.0        1.0    
   1.00023    2.00001    1.00002
   1.00046    2.00002    1.00004
   1.00068    2.00003    1.00006
   1.00091    2.00004    1.00008
   1.00114    2.00005    1.0001 
   1.00137    2.00006    1.00012
   1.00159    2.00007    1.00014
   1.00182    2.00008    1.00016
   1.00205    2.00009    1.00018
   1.00228    2.0001     1.0002 
   1.0025     2.00011    1.00022
   1.00273    2.00012    1.00024
   ⋮                            
 NaN        NaN        NaN      
 NaN        NaN        NaN      
 NaN        NaN        NaN      
 NaN        NaN        NaN      
 NaN        NaN        NaN      
 NaN        NaN        NaN      
 NaN        NaN        NaN      
 NaN        NaN        NaN      
 NaN        NaN        NaN      
 NaN        NaN        NaN      
 NaN        NaN        NaN      
 NaN        NaN        NaN      

Runge-Kutta de orden 2

Primero para el intervalo [1,1.8]:

In [15]:
tRK2_1 = collect(1:0.0001:1.8)
solucion_RK2_1 = RK_2(edo_y,cond_inic_y,tRK2_1)
Out[15]:
8001×3 Array{Float64,2}:
   1.0       2.0      1.0    
   1.00228   2.0001   1.0002 
   1.00455   2.0002   1.0004 
   1.00683   2.0003   1.0006 
   1.0091    2.0004   1.0008 
   1.01138   2.0005   1.001  
   1.01365   2.0006   1.0012 
   1.01593   2.00071  1.0014 
   1.01821   2.00081  1.0016 
   1.02048   2.00091  1.0018 
   1.02276   2.00101  1.002  
   1.02503   2.00111  1.0022 
   1.02731   2.00122  1.0024 
   ⋮                         
 448.694    41.1314   7.38565
 449.438    41.1763   7.38976
 450.183    41.2212   7.39388
 450.93     41.2663   7.39801
 451.679    41.3114   7.40213
 452.43     41.3566   7.40627
 453.182    41.4019   7.41041
 453.935    41.4473   7.41455
 454.691    41.4927   7.41869
 455.448    41.5382   7.42285
 456.207    41.5838   7.427  
 456.967    41.6295   7.43116

Y ahora para el intervalo [1,2.1]:

In [16]:
tRK2_2 = collect(1:0.0001:2.1)
solucion_RK2_2= RK_2(edo_y,cond_inic_y,tRK2_2)
DomainError with -Inf:
cos(x) is only defined for finite x.

Stacktrace:
 [1] cos_domain_error(::Float64) at .\special\trig.jl:98
 [2] cos(::Float64) at .\special\trig.jl:109
 [3] edo_y(::Array{Float64,1}, ::Float64) at .\In[9]:3
 [4] RK_2(::typeof(edo_y), ::Array{Int64,1}, ::Array{Float64,1}) at .\In[5]:7
 [5] top-level scope at In[16]:2

Runge-Kutta de orden 4

Primero para el intervalo [1,1.8]:

In [17]:
tRK4_1 = collect(1:0.0001:1.8)
solucion_RK4_1 = RK_4(edo_y,cond_inic_y,tRK4_1)
Out[17]:
8001×3 Array{Float64,2}:
   1.0       2.0      1.0    
   1.00228   2.0001   1.0002 
   1.00455   2.0002   1.0004 
   1.00683   2.0003   1.0006 
   1.0091    2.0004   1.0008 
   1.01138   2.0005   1.001  
   1.01365   2.0006   1.0012 
   1.01593   2.00071  1.0014 
   1.01821   2.00081  1.0016 
   1.02048   2.00091  1.0018 
   1.02276   2.00101  1.002  
   1.02503   2.00111  1.0022 
   1.02731   2.00122  1.0024 
   ⋮                         
 448.694    41.1314   7.38565
 449.438    41.1763   7.38976
 450.184    41.2213   7.39388
 450.931    41.2663   7.39801
 451.68     41.3115   7.40214
 452.43     41.3567   7.40627
 453.182    41.402    7.41041
 453.936    41.4473   7.41455
 454.691    41.4927   7.4187 
 455.449    41.5382   7.42285
 456.207    41.5838   7.42701
 456.968    41.6295   7.43117

Y ahora para el intervalo [1,2.1]:

In [18]:
tRK4_2 = collect(1:0.0001:2.1)
solucion_RK4_2 = RK_4(edo_y,cond_inic_y,tRK4_2)
DomainError with -Inf:
cos(x) is only defined for finite x.

Stacktrace:
 [1] cos_domain_error(::Float64) at .\special\trig.jl:98
 [2] cos(::Float64) at .\special\trig.jl:109
 [3] edo_y(::Array{Float64,1}, ::Float64) at .\In[9]:3
 [4] RK_4(::typeof(edo_y), ::Array{Int64,1}, ::Array{Float64,1}) at .\In[6]:6
 [5] top-level scope at In[18]:2

En todas las soluciones que obtuvimos, observamos que todos los métodos tienen problemas a partir de aproximadamente 1.8, ya que ahí, se presentan indeterminaciones.

Finalmente, vamos a graficar la solución con el método de Rugen-Kutta de orden 4 y con el método de Euler con Jacobiano numérico para poder visualizar su comportamiento cerca de la zona de indeterminación que detectamos arriba, así como también, para comprobar que el algoritmo desarrollado en el problema 1 funciona correctamente:

In [19]:
plot(title="Solución en el intervalo [1, 1.8]",xlabel="t") 
plot!(t_euler_n1,solucion_euler_n1[:,1],label="g_1(t) con Euler")
scatter!(tRK4_1,solucion_RK4_1[:,1],m=:cross,label="g_1(t) con RK4")
plot!(t_euler_n1,solucion_euler_n1[:,2],label="g_2(t) con Euler")
scatter!(tRK4_1,solucion_RK4_1[:,2],m=:cross,label="g_2(t) con RK4")
plot!(t_euler_n1,solucion_euler_n1[:,3],label="g_3(t) con Euler")
scatter!(tRK4_1,solucion_RK4_1[:,3],m=:cross,label="g_3(t) con RK4")
Out[19]:

Primero, notamos que las soluciones con ambos métodos son similares, por lo que, podemos concluir que el algoritmo del problema 1 trabaja adecuadamente. Por otra parte, observando la gráfica vemos que efectivamente tenemos una variable que diverge y por tanto, no podemos obtener una solución acertada en dicha zona con ninguno de los métodos. De hecho, para asegurarnos aún más, podemos tomar nuestro intervalo de tiempo como la unión de dos intervalos, en donde el segundo intervalo va a ir desde 1.7 hasta 2.1 pero con muchos más puntos que en el primer intervalo para tratar de visualizar mejor el comportamiento justo en la zona donde habíamos visto que tenemos problemas:

In [20]:
t_euler_n21 = collect(1.7:0.08:2.1) # t entre 1.7 y 2.1
time=vcat(t_euler_n1,t_euler_n21)
solucion_analitico_21 = Euler_mejorado_J_t_p(edo_y_1,cond_inic_y,time,jacob_edo_y)
#solucion_numerico_21=Euler_mejorado_Jacobnum(edo_y,edo_y_1,cond_inic_y,time)
plot(title="Solución en el intervalo [1, 2.1]",xlabel="t") 
plot!(time,solucion_analitico_21[:,1],label="g_1(t)")
plot!(time,solucion_analitico_21[:,2],label="g_2(t)")
plot!(time,solucion_analitico_21[:,3],label="g_3(t)")
Out[20]:

Efectivamente hay problemas a partir de 1.75 aprox, y de hecho ni siquiera logra graficarlo bien.

3. Modifique las funciones integradoras que programamos en clase (Euler mejorado, Runge-Kutta de $2^°$ orden y Runge-Kutta de $4^°$ orden) para que se puedan pasar argumentos opcionales al sistema de ecuaciones diferenciales de tal manera que se puedan modificar parámetros al momento de hacer la integración. Una vez hecho esto, obtenga la solución de la siguiente ecuación diferencial:

$$ \ddot{x} + \frac{1}{10} \dot{x} + 4 sin x = \frac{1}{2} sin (\omega t) $$

Con condiciones iniciales $\dot{x}(0) = 1$ y $x(0) = 0$ y para $\omega = {1, 1.1, 1.2, 1.3, ..., 2.9, 3}$. Realice la gráfica de algunas soluciones que representen el comportamiento de la solución en función de $\omega$ y explique su resultado. La solución debe tener una precisión de por lo menos $10^{−5}$.

Euler mejorado

In [21]:
function Euler_mejorado_J_t_p_omega(edo,p_ini,t,jacob,ω) 
    #anexamos omega, que va a ser el parámetro que vamos a poder modificar
    # J y g con dependencia temporal
    sol = zeros(length(t),length(p_ini))
    sol[1,:] = p_ini
    δ = t[2]-t[1]
    for i in 1:(length(t)-1)
        eval_edo = edo(sol[i,:],t[i],ω) #Añadimos ω porque ahora la edo depende de ella
        if length(eval_edo) == length(p_ini)
            sol[i+1,:] .= sol[i,:] .+ δ .*eval_edo
        else
            sol[i+1,:] .= sol[i,:] .+ δ .*eval_edo[1:(end-1)]
        end
        sol[i+1,:] .+= 0.5*(δ^2) .*(jacob(sol[i,:],t[i],ω) *eval_edo )
    end
    return sol
end
Out[21]:
Euler_mejorado_J_t_p_omega (generic function with 1 method)

Runge-Kutta de $2^°$ orden

In [22]:
function RK_2_omega(edo,x_ini,t,ω)
    #anexamos ω, que va a ser el parámetro que vamos a poder modificar
    sol = zeros( length(t) , length(x_ini) )
    sol[1,:] .= x_ini
    δ = t[2]-t[1]
    for i in 1:length(t)-1
        k1 = sol[i,:] .+ 0.5*δ .*edo(sol[i,:],t[i],ω)
        sol[i+1,:] .= sol[i,:] .+ δ*edo(k1,t[i]+0.5*δ,ω)
    end
    return sol
end
Out[22]:
RK_2_omega (generic function with 1 method)

Runge-Kutta de $4^°$ orden

In [23]:
function RK_4_omega(edo,x_ini,t,ω)
    #anexamos ω, que va a ser el parámetro que vamos a poder modificar
    sol = zeros( length(t) , length(x_ini) )
    sol[1,:] .= x_ini
    δ = t[2]-t[1]
    for i in 1:length(t)-1
        k1 = edo(sol[i,:],t[i],ω)
        k2 = edo(sol[i,:] .+ 0.5*δ.*k1 , t[i] + 0.5*δ,ω)
        k3 = edo(sol[i,:] .+ 0.5*δ.*k2 , t[i] + 0.5*δ,ω)
        k4 = edo(sol[i,:] .+ δ.*k2 , t[i] + δ,ω)
        sol[i+1,:] .= sol[i,:] .+ (δ/6.0).*(k1 .+ 2.0.*k2 .+ 2.0.*k3 .+ k4)
    end
    return sol
end
Out[23]:
RK_4_omega (generic function with 1 method)

EDO y Jacobiano

Vamos a analizar ahora la ecuación diferencial que tenemos:

Tenemos:

$$ \ddot{x} + \frac{1}{10} \dot{x} + 4 sin x = \frac{1}{2} sin (\omega t) $$

Con condiciones iniciales $\dot{x}(0) = 1$ y $x(0) = 0$ y para $\omega = {1, 1.1, 1.2, 1.3, ..., 2.9, 3}$.

Siguiendo el análisi realizado en la ecuación del problema 2, podemos ver a $x$ como la posición:

$$ \Rightarrow \dot{x}= v \quad (\text{velocidad}) $$$$ \Rightarrow \ddot{x}=\dot{v} \quad (\text{aceleración}) $$

Por lo que, la ecuación diferencial se transforma en:

$$ \dot{v} + \frac{1}{10} v + 4 sin (x) = \frac{1}{2} sin (\omega t) $$

Y tenemos las siguientes 2 EDO's de 1er orden:

$$ \dot{v} = \frac{1}{2} sin (\omega t) - \frac{1}{10} v - 4 sin (x) $$$$ \dot{x} = v $$

Y utilizando las condiciones iniciales dadas sabemos que: $v(0)=1$, $x(0)=0$

Por lo tanto, podemos definir las funciones g's como:

$$ g_1=\frac{1}{2} sin (\omega t) - \frac{1}{10} v - 4 sin (x) $$$$ g_2 = v $$

A continuación, procedemos a programar la ecuación diferencial:

In [24]:
function edo_x(variables,t,ω) 
    #variables=(v,x)
    g1=(0.5*sin(ω*t))-((1/10)*variables[1]) - (4*sin(variables[2]))
    g2=variables[1]
    return [g1,g2,1]
end
Out[24]:
edo_x (generic function with 1 method)

Definimos las condiciones iniciales y el intervalo de tiempo que ocuparemos:

In [25]:
cond_inic_x=[1,0]
tiempo=collect(0:0.1:10*π);

Y calculamos su Jacobiano de manera analítica:

$$\mathbb{J}^* = \begin{bmatrix} \frac{\partial g_1}{\partial v} & \frac{\partial g_1}{\partial x} & \frac{\partial g_1}{\partial t} \\ \frac{\partial g_2}{\partial v} & \frac{\partial g_2}{\partial x} & \frac{\partial g_2}{\partial t} \\ \end{bmatrix}$$$$ \Rightarrow \mathbb{J}^* = \begin{bmatrix} -\frac{1}{10} & -4 cos(x) & \frac{1}{2}\omega cos(\omega t) \\ 1 & 0 & 0 \end{bmatrix}$$

A continuación, procedemos a programar este Jacobiano:

In [26]:
function jacob_edo_x(variables,t,ω)
    Jacob = zeros(2,3) #definimos una matriz de ceros de 2x3 para poder llenarla
    Jacob[1,1] = -1/10
    Jacob[1,2] = -4*cos(variables[2])
    Jacob[1,3] = 0.5*ω*cos(ω*t)
    Jacob[2,1] = 1
    return Jacob
end
Out[26]:
jacob_edo_x (generic function with 1 method)

Y por último, recordamos que para los métodos de Runge-Kutta debemos volver a definir la EDO pero sin agregarle el 1 extra al arreglo que regresa la función.

In [27]:
function edo_x_RK(variables,t,ω) 
    #variables=(v,x)
    g1=(0.5*sin(ω*t))-((1/10)*variables[1]) - (4*sin(variables[2]))
    g2=variables[1]
    return [g1,g2]
end
Out[27]:
edo_x_RK (generic function with 1 method)

Con esto, ya tenemos todo lo necesario para obtener la solución con los 3 métodos. Ahora, el valor de $\omega$ no influye en la obtención de las cifras significativas ya que se trata de un parámetro constante, por tanto, basta con ver que se obtienen las 5 cifras significativas para una ω en particular. Vamos a obtener las soluciones para $\omega=1.1$:

In [28]:
solucion_euler_ω = Euler_mejorado_J_t_p_omega(edo_x,cond_inic_x,tiempo,jacob_edo_x,1.1) #solucion con el metodo de Euler
Out[28]:
315×2 Array{Float64,2}:
  1.0        0.0       
  0.9728     0.0995    
  0.91242    0.194581  
  0.822022   0.282046  
  0.706003   0.359081  
  0.569575   0.423364  
  0.41836    0.473127  
  0.258067   0.507173  
  0.0942823  0.524877  
 -0.0676317  0.526163  
 -0.222571   0.511479  
 -0.365698   0.481772  
 -0.492474   0.438457  
  ⋮                    
 -0.174011   0.103502  
 -0.164017   0.0864759 
 -0.149969   0.0706764 
 -0.132995   0.0564559 
 -0.114295   0.0440488 
 -0.0950871  0.0335673 
 -0.0765538  0.025002  
 -0.0597875  0.0182288 
 -0.0457431  0.0130197 
 -0.0351973  0.00905929
 -0.0287157  0.0059642 
 -0.0266302  0.00330564
In [29]:
solucion_RK2_ω = RK_2_omega(edo_x_RK,cond_inic_x,tiempo,1.1) #solucion con el metodo de Runge-Kutta de orden 2
Out[29]:
315×2 Array{Float64,2}:
  1.0        0.0       
  0.972807   0.0995    
  0.912472   0.194582  
  0.822141   0.282051  
  0.706189   0.359098  
  0.569806   0.4234    
  0.418599   0.473185  
  0.258272   0.507254  
  0.0944166  0.524977  
 -0.0675923  0.526274  
 -0.222634   0.511593  
 -0.365855   0.481877  
 -0.492704   0.438545  
  ⋮                    
 -0.174312   0.103369  
 -0.164325   0.0863153 
 -0.150271   0.0704883 
 -0.133276   0.0562415 
 -0.11454    0.0438107 
 -0.0952822  0.0333096 
 -0.0766867  0.0247301 
 -0.0598476  0.0179491 
 -0.0457222  0.0127396 
 -0.03509    0.00878687
 -0.0285196  0.00570791
 -0.0263462  0.00307398
In [30]:
solucion_RK4_ω = RK_4_omega(edo_x_RK,cond_inic_x,tiempo,1.1) #solucion con el metodo de Runge-Kutta de orden 4
Out[30]:
315×2 Array{Float64,2}:
  1.0         0.0       
  0.972938    0.0989269 
  0.913185    0.193488  
  0.823816    0.280556  
  0.709094    0.357377  
  0.574084    0.421663  
  0.424269    0.47166   
  0.265241    0.506166  
  0.1025      0.524538  
 -0.0586631   0.526672  
 -0.213199    0.51298   
 -0.356324    0.484368  
 -0.483551    0.442204  
  ⋮                     
 -0.210295    0.133614  
 -0.211206    0.112486  
 -0.206171    0.091571  
 -0.195953    0.0714265 
 -0.181485    0.0525254 
 -0.163819    0.0352409 
 -0.144083    0.0198365 
 -0.123428    0.00646139
 -0.102974   -0.00484909
 -0.0837655  -0.0141683 
 -0.0667243  -0.0216682 
 -0.0526135  -0.0276054 

Observamos que para los 3 métodos se cumple que se tienen cinco cifras significativas como minímo. Por lo tanto, procedemos ahora sí a realizar las gráficas para presentar las soluciones con distintas ω:

In [31]:
ωs=collect(1:0.1:3)
plot(title="Solución para v(t) con distintos valores de ω (Euler)",xlabel="t",ylabel="v(t)")
for ω_i in ωs
    solucion = Euler_mejorado_J_t_p_omega(edo_x,cond_inic_x,tiempo,jacob_edo_x,ω_i)
    plot!(tiempo,solucion[:,1],label ="ω=$(ω_i)", lw=1.0)
end
plot!()
Out[31]:
In [32]:
plot(title="Solución para x(t) con distintos valores de ω (Euler)",xlabel="t",ylabel="x(t)")
for ω_i in ωs
    solucion = Euler_mejorado_J_t_p_omega(edo_x,cond_inic_x,tiempo,jacob_edo_x,ω_i)
    plot!(tiempo,solucion[:,2],label ="ω=$(ω_i)", lw=1.0)
end
plot!()
Out[32]:
In [33]:
plot(title="Solución para v(t) con distintos valores de ω (Runge-Kutta 2)",xlabel="t",ylabel="v(t)")
for ω_i in ωs
    solucion = RK_2_omega(edo_x_RK,cond_inic_x,tiempo,ω_i)
    plot!(tiempo,solucion[:,1],label ="ω=$(ω_i)", lw=1.0)
end
plot!()
Out[33]:
In [34]:
plot(title="Solución para x(t) con distintos valores de ω (Runge-Kutta2)",xlabel="t",ylabel="x(t)")
for ω_i in ωs
    solucion = RK_2_omega(edo_x_RK,cond_inic_x,tiempo,ω_i)
    plot!(tiempo,solucion[:,2],label ="ω=$(ω_i)", lw=1.0)
end
plot!()
Out[34]:
In [35]:
plot(title="Solución para v(t) con distintos valores de ω (Runge-Kutta 4)",xlabel="t",ylabel="v(t)")
for ω_i in ωs
    solucion = RK_4_omega(edo_x_RK,cond_inic_x,tiempo,ω_i)
    plot!(tiempo,solucion[:,1],label ="ω=$(ω_i)", lw=1.0)
end
plot!()
Out[35]:
In [36]:
plot(title="Solución para x(t) con distintos valores de ω (Runge-Kutta4)",xlabel="t",ylabel="x(t)")
for ω_i in ωs
    solucion = RK_4_omega(edo_x_RK,cond_inic_x,tiempo,ω_i)
    plot!(tiempo,solucion[:,2],label ="ω=$(ω_i)", lw=1.0)
end
plot!()
Out[36]:

Podemos notar que tanto para $x(t)$, como $v(t)$, tenemos un fenómeno de resonancia. Asimismo, conforme el sistema va avanzando en el tiempo, la amplitud de las soluciones va creciendo, haciendo que mientras mayor sea el valor de $\omega$ también sea mayor el valor de la amplitud.

4. Considere el siguiente sistema de ecuaciones diferenciales:

$$ \frac{dx}{dt} = 10 (y-x) $$

$$ \frac{dy}{dt} = x (28-z)-y $$

$$ \frac{dz}{dt} = xy - \frac{8}{3} z $$

Realice una integración numérica del sistema de ecuaciones con el método de Runge-Kutta de $4^°$ y realice los siguientes ejercicios:

a) Elija una condición inicial aleatoria en el $x, y, z \in [−5, 5]$ y realice la gráfica de la solución en el plano xy, xz y yz y comente sus observaciones.

b) Ahora tome 3 soluciones parecidas. Para ello proponga 3 condiciones iniciales tales que en la coordenada en $y$ las condiciones iniciales difieran entre si en $10^{−2}$. Explique lo que observa.

Nota: De preferencia realice integraciones con tiempos relativamente largos en los que pueda apreciar diferencias en el comportamiento de las soluciones. Además asegure que cada una de las soluciones que encuentre sea convergente en por lo menos 5 cifras ($10^{−5}$), es decir encuentre un paso adecuado para realizar la integración.

Vamos a definir primero nuestro sistema de ecuaciones:

In [37]:
function edo_xyz(variables,t)
    #variables(x,y,z)
    dx=10*(variables[2]-variables[1])
    dy=variables[1]*(28-variables[3])-variables[2]
    dz=variables[1]*variables[2]-((8/3)*variables[3])
    return[dx,dy,dz]
end
Out[37]:
edo_xyz (generic function with 1 method)

a) Utilizaremos la función rand para elegir una condición inicial aleatoria:

In [38]:
cond_inic_xyz=[-5rand(),5rand(),5rand()];

Gráficas de la solución:

In [39]:
plot(title="Solución de la EDO (Plano xy)",xlabel="x",ylabel="y")
tiempo = collect(0:0.01:60π)
sol_num_rk4 = RK_4(edo_xyz,cond_inic_xyz,tiempo)
plot!(sol_num_rk4[:,1],sol_num_rk4[:,2])
plot!(legend=false)
Out[39]:
In [40]:
plot(title="Solución de la EDO (Plano xz)",xlabel="x",ylabel="z")
tiempo = collect(0:0.01:60π)
sol_num_rk4 = RK_4(edo_xyz,cond_inic_xyz,tiempo)
plot!(sol_num_rk4[:,1],sol_num_rk4[:,3])
plot!(legend=false)
Out[40]:
In [41]:
plot(title="Solución de la EDO (Plano yz)",xlabel="y",ylabel="z")
tiempo = collect(0:0.01:60π)
sol_num_rk4 = RK_4(edo_xyz,cond_inic_xyz,tiempo)
plot!(sol_num_rk4[:,2],sol_num_rk4[:,3])
plot!(legend=false)
Out[41]:

En las 3 gráficas obtenemos la misma forma pero vista desde distintos ángulos. Asimismo, podemos ver que conforme el tiempo va creciendo, aparecen más líneas que siguen la misma trayectoria.

Finalmente, vamos a comprobar que la solución tiene al menos 5 cifras significativas:

In [42]:
sol_num_rk4 = RK_4(edo_xyz,cond_inic_xyz,tiempo)
Out[42]:
18850×3 Array{Float64,2}:
 -4.62825   3.61768     1.69786
 -3.90042   2.46561     1.52564
 -3.34279   1.48806     1.4151 
 -2.92448   0.64572     1.34499
 -2.62101  -0.0946211   1.3021 
 -2.41313  -0.760892    1.27846
 -2.28586  -1.37663     1.26955
 -2.22769  -1.96185     1.27326
 -2.23     -2.53374     1.28912
 -2.28651  -3.10736     1.31797
 -2.39286  -3.69613     1.36173
 -2.54637  -4.31226     1.42338
 -2.74569  -4.96707     1.50696
  ⋮                            
 -2.44071  -0.174169   24.4034 
 -2.22927  -0.263235   23.7663 
 -2.04667  -0.357105   23.1474 
 -1.89055  -0.454344   22.5462 
 -1.75867  -0.554001   21.9619 
 -1.64892  -0.655534   21.3942 
 -1.55935  -0.758731   20.8424 
 -1.48821  -0.863657   20.3061 
 -1.43392  -0.970601   19.785  
 -1.39509  -1.08004    19.2786 
 -1.37052  -1.19258    18.7868 
 -1.35918  -1.309      18.3093 

Efectivamente obtuvimos una solución con al menos 5 cifras significativas.

b) Ahora, para tomar 3 condiciones iniciales aleatorias y parecidas realizaremos lo siguiente:

In [43]:
#Partimos de la condicion inicial que ya teníamos:
cond_inic_1=[cond_inic_xyz[1],cond_inic_xyz[2],cond_inic_xyz[3]]
#Y para tener otras 2 condiciones iniciales cercanas, le vamos a sumar y restar 0.01 para que la diferencia entre ellas sea de 10^{-2}
cond_inic_2=[cond_inic_xyz[1],cond_inic_xyz[2]+0.01,cond_inic_xyz[3]]
cond_inic_3=[cond_inic_xyz[1],cond_inic_xyz[2]-0.01,cond_inic_xyz[3]]
Out[43]:
3-element Array{Float64,1}:
 -4.628248920593543 
  3.6076769048224246
  1.6978628734426915
In [44]:
plot(title="3 Soluciones parecidas de la EDO (Plano xy)",xlabel="x",ylabel="y")
tiempo = collect(0:0.01:60π)
solnum_rk4_1 = RK_4(edo_xyz,cond_inic_1,tiempo)
solnum_rk4_2 = RK_4(edo_xyz,cond_inic_2,tiempo)
solnum_rk4_3 = RK_4(edo_xyz,cond_inic_3,tiempo)
plot!(solnum_rk4_1[:,1],solnum_rk4_1[:,2], label="Cond. inicial 1")
plot!(solnum_rk4_2[:,1],solnum_rk4_2[:,2], label="Cond. inicial 2")
plot!(solnum_rk4_3[:,1],solnum_rk4_3[:,2], label="Cond. inicial 3")
plot!(legend=true, grid=true)
Out[44]:
In [45]:
plot(title="3 Soluciones parecidas de la EDO (Plano xz)",xlabel="x",ylabel="y")
tiempo = collect(0:0.01:60π)
solnum_rk4_1 = RK_4(edo_xyz,cond_inic_1,tiempo)
solnum_rk4_2 = RK_4(edo_xyz,cond_inic_2,tiempo)
solnum_rk4_3 = RK_4(edo_xyz,cond_inic_3,tiempo)
plot!(solnum_rk4_1[:,1],solnum_rk4_1[:,3], label="Cond. inicial 1")
plot!(solnum_rk4_2[:,1],solnum_rk4_2[:,3], label="Cond. inicial 2")
plot!(solnum_rk4_3[:,1],solnum_rk4_3[:,3], label="Cond. inicial 3")
plot!(legend=true, grid=true)
Out[45]:
In [46]:
plot(title="3 Soluciones parecidas de la EDO (Plano yz)",xlabel="x",ylabel="y")
tiempo = collect(0:0.01:60π)
solnum_rk4_1 = RK_4(edo_xyz,cond_inic_1,tiempo)
solnum_rk4_2 = RK_4(edo_xyz,cond_inic_2,tiempo)
solnum_rk4_3 = RK_4(edo_xyz,cond_inic_3,tiempo)
plot!(solnum_rk4_1[:,2],solnum_rk4_1[:,3], label="Cond. inicial 1")
plot!(solnum_rk4_2[:,2],solnum_rk4_2[:,3], label="Cond. inicial 2")
plot!(solnum_rk4_3[:,2],solnum_rk4_3[:,3], label="Cond. inicial 3")
plot!(legend=true, grid=true)
Out[46]:

Observando las gráficas anteriores, podemos notar que todas las soluciones tienen la misma forma sin importar la condición inicial, sin embargo si se puede apreciar un ligero desfase entre ellas (ya que se pueden distinguir los 3 colores), por lo que no podemos decir que son iguales, lo cual tiene sentido ya que las soluciones dependen justo de las condiciones iniciales que se tengan.

Finalmente, vamos a comprobar que para cada una de las soluciones tengamos al menos 5 cifras significativas:

In [47]:
solnum_rk4_1 = RK_4(edo_xyz,cond_inic_1,tiempo)
Out[47]:
18850×3 Array{Float64,2}:
 -4.62825   3.61768     1.69786
 -3.90042   2.46561     1.52564
 -3.34279   1.48806     1.4151 
 -2.92448   0.64572     1.34499
 -2.62101  -0.0946211   1.3021 
 -2.41313  -0.760892    1.27846
 -2.28586  -1.37663     1.26955
 -2.22769  -1.96185     1.27326
 -2.23     -2.53374     1.28912
 -2.28651  -3.10736     1.31797
 -2.39286  -3.69613     1.36173
 -2.54637  -4.31226     1.42338
 -2.74569  -4.96707     1.50696
  ⋮                            
 -2.44071  -0.174169   24.4034 
 -2.22927  -0.263235   23.7663 
 -2.04667  -0.357105   23.1474 
 -1.89055  -0.454344   22.5462 
 -1.75867  -0.554001   21.9619 
 -1.64892  -0.655534   21.3942 
 -1.55935  -0.758731   20.8424 
 -1.48821  -0.863657   20.3061 
 -1.43392  -0.970601   19.785  
 -1.39509  -1.08004    19.2786 
 -1.37052  -1.19258    18.7868 
 -1.35918  -1.309      18.3093 
In [48]:
solnum_rk4_2 = RK_4(edo_xyz,cond_inic_2,tiempo)
Out[48]:
18850×3 Array{Float64,2}:
 -4.62825    3.62768     1.69786
 -3.89947    2.47563     1.52524
 -3.34096    1.49833     1.41437
 -2.92183    0.65645     1.34398
 -2.61756   -0.0832194   1.30083
 -2.40888   -0.748617    1.2769 
 -2.28079   -1.36328     1.2677 
 -2.22178   -1.94722     1.27105
 -2.22319   -2.51762     1.28649
 -2.27873   -3.08953     1.31483
 -2.38404   -3.67635     1.35797
 -2.5364    -4.29029     1.41884
 -2.73446   -4.94266     1.50145
  ⋮                             
  4.92707   -1.06474    30.4445 
  4.35247   -1.14752    29.5926 
  3.82698   -1.1843     28.7668 
  3.34978   -1.18583    27.9678 
  2.91917   -1.16119    27.1956 
  2.53284   -1.1179     26.4493 
  2.18803   -1.0621     25.7279 
  1.88175   -0.998674   25.0302 
  1.61086   -0.931474   24.355  
  1.3722    -0.863456   23.7009 
  1.16267   -0.79684    23.0669 
  0.979288  -0.733244   22.4518 
In [49]:
solnum_rk4_3 = RK_4(edo_xyz,cond_inic_3,tiempo)
Out[49]:
18850×3 Array{Float64,2}:
 -4.62825   3.60768   1.69786
 -3.90137   2.4556    1.52605
 -3.34461   1.47779   1.41583
 -2.92713   0.63499   1.346  
 -2.62446  -0.106023  1.30338
 -2.41738  -0.773167  1.28001
 -2.29092  -1.38998   1.27141
 -2.2336   -1.97648   1.27547
 -2.23681  -2.54986   1.29175
 -2.29428  -3.12519   1.32111
 -2.40169  -3.7159    1.36551
 -2.55634  -4.33422   1.42794
 -2.75692  -4.99148   1.51249
  ⋮                          
  1.57608   2.92629   9.44347
  1.71773   3.20273   9.24475
  1.87355   3.50744   9.06091
  2.04501   3.8432    8.8935 
  2.23374   4.21306   8.74448
  2.44148   4.62028   8.61621
  2.67015   5.06833   8.51161
  2.92183   5.56088   8.43421
  3.19873   6.10173   8.3883 
  3.50324   6.69475   8.37906
  3.83788   7.34379   8.41273
  4.20531   8.05248   8.49676

Efectivamente, para los 3 casos tenemos una solución con al menos 5 cifras significativas.